- PREPRINT -

How to build a real-time prediction model for influenza in Germany

Paul Schneider, Maastricht University, Netherlands Institute of Health Service Research
David Barnett, Maastricht University
John Paget, Netherlands Institute of Health Service Research
Peter Spreeuwenberg, Netherlands Institute of Health Service Research
Christel van Gool, Maastricht University

Contact: schneider.paulpeter@gmail.com

---

Summary

This is a fully reproducable R-tutorial on how to build a (near) real-time prediction model for influenza in Germany, based on search query data from Google Trends and page view statistics from Wikipedia. All code and data is openly available. Instructions include guidance on how to identify and download relevant data, format and pre-process the retrieved information, and build and evaluate prediction models, using multiple statistical modelling techniques (a.k.a. machine learning). In addition, we will provide a short introduction to the emerging field of ‘Digital Epidemiology’ and discuss some of its critical methodological issues.

TL;DR: Code only version

Methods are inspired by the caret R-package and the book Applied Predictive Modelling, both by Max Kuhn.

****

Introduction

Digital Epidemiology

In the past, people could only went to see their doctor when they were sick and wanted to get information, guidance and treatment. Many national disease surveillance systems are therefore based on doctors as sentinels: health care providers and medical laboratories are being asked to report the number of cases for a certain disease to the surveillance organization. In Germany, for example, laboratories are required to report any proof of influenza viruses to the local health authority, which in turn reports it to the Robert Koch Institute. Nowadays, however, many people first turn to the internet. They seek health information by googling for their sysmptoms, look up diseases or treatments on Wikipedia or connect with others through social media. For Epidemiologists, this offers the opportunity to track these digital traces of people. This has been described as ‘Digital Epidemilogy’: “Digital data sources, when harnessed appropriately, can provide local and timely information about disease and health dynamics in populations around the world”.

Milestones of digital epidemiology

• In 2008, Google Flu Trend was launched and accompanied by a broad media coverage. ..
• In 2014, McIver and Brownstein claimed that Wikipedia page view data …
• Twitter
• Health Map
• Influenza mao netherlands

Why Influenza?

In principle, the methods of surveillance, we use in this example, can be applied to any disease. However, some of the features of influenza make it easier to detect relevant predictors and monitor its activity. First, influenza has a strong seasonality (in most countries). Every year there is a start and an end, and during the summer months there are literally no cases of influenza. This strong annual pattern can be detected easily. In diseases, in which the incidence is more stable, there might be not enough variation to distinguish signal and noise. Second, during the flu season, between 5% and 20% of the population get sick. The online traces of rare diseases (say, achalasia) probably get lost in noise. And finally, influenza has distinct symptoms that most people will recognize as flu-symptoms, this means the symptoms are (somewhat) specific to influenza. In contrast, diseases with only subtle or very non-specific symptoms (say, fatigue) are certainly more problematic to track.

Nevertheless, there is no reason to assume that other diseases can not be surveyed using similar methods. A few papers have already shown promising results in …Ebola, TUberculosis… [ADD LINK][LINK]. We would indeed be very interested to see how predictive models perform in further diseases, like allergic rhinitis, depression or even sunburns. If you have good data on this or know where to get it, we invite you to follow this tutorial, build your own model(s) and share your results with us and the public.

Issues

The huge amount of available data and the novelity of methods are a challenge.
Overly optimistic because of

• Variety of methods and data sources.
• Inconsistent reporting of methods and results, the ai should be reproducability. In the area of digital epidemiology, this might imply open sharing of code and data- otherwise reproduciability is rarely possible.
• The failure to strictly separate trainig and test data sets
  • Espeically in split-cross validation, leaking info from the future to the past
• reporting of correlation coeficients, adjusted or not, or rmse, both?!
• Research area mostly in the US (cite global paper), unclear whether methods are reproducible in countries in which digitalization is not that advanced.

Research within digital epi must adhere to the rigour and caution other epidemiological project are also held accountable to.

---

Tutorial

Overview

We will build a model for predicting the influenza incidence (lab-confirmed cases) in Germany in near real-time (‘Nowcasting’), by using data from various sources. We will start with getting information on influenza incidence in Germany, since the availability of the outcome data determines the scope and time horizon of the model we build. We then download data on search queries and page view statistics from Google Correlate, Google Trends and Wikipedia, respectively. The data have to be formatted, merged and pre-processed, before entering into the model. Subsequently, we will use various so called machine learning methods to select relevant predictors and build prediction models. Finally, we will compare the different models in terms of predictive accuracy and select the best model to predict this week’s influenza activity.

We would like to point out that when we speak of ‘predictions’, we use the term in its statistical sense. We do not mean ’forecasting, as in forecasting future influenza cases. We are doing ‘Nowcasting’: We observe what people were doing on the internet until today and infer from this information how many influenza cases there are about now

Data sources

Outome data is taken from the Robert Koch Institute, predictor data from Google Correlate, Google Trends and Wikipedia via it’s API and Wikishark. All analyses are performed using the statistical software R/R Studio, and plenty of its packages, particularly the caret package. Everything we use is available online for free. The data used in this example is also available on github. If you want to build your own model, you have to have a Google account to upload outcome data to Google Correlate.

Transferability

If you want to transfer this approach to another setting (another country or another disease), the only thing you need is a set of outcome data, i.e. good-quality data of the actual weekly incidence of the disease, over a sufficient time span (1 year +), in a specific country. The rest should work (more or less) automatically, unless you wish to add more data sources as predictors, which may require some more editing. Depending on the country in which you want to predict disease activity, the availability and quality of Google and Wikipedia data may differ substantially. For Wikipedia data, page view statistics can only be distinguished by language, but not by country. Language can be a good proxy in some cases (e.g. Japanese, Korean, Italian, etc.), but a pretty bad one in others (e.g. French, English, Spanish). You can look up how much of a country’s Wikipedia traffic is in a specfific language, and how much Wikipedia traffic within a specific language is from a specfific country. Moreover, with weekly data, you should be aware that there are divergent opinions about what is the best starting day for a week (Saturday? Sunday? Monday?).

---

Getting started

Before building the model, we need to install and load the required R-packages (yes, there are a lot)…

# Install and load all required packages
# This may take a while...
required_packages<-c("knitr","RCurl","ISOweek","jsonlite","ggplot2","prophet","dplyr","gtrendsR","wikipediatrend","pageviews","caret","imputeTS","gridExtra","cowplot","corrplot","doParallel","glmnet", "Cubist","pls","devtools","plotmo","plotrix","TeachingDemos","earth","kernlab","rpart","party","grid","zoo","mvtnorm","modeltools","stats4","strucchange","sandwich","elasticnet","lars","randomForest","ipred","e1071") # 
pft_packages <- function(package){
    for(i in 1:length(package)){
      if(eval(parse(text=paste("require(",package[i],")")))==0) {
        install.packages(package)}}
    return (eval(parse(text=paste("require(",package,")"))))}
pft_packages(required_packages)
#  gtrendR has to be the developer version. If you don't already have it, install it throughgithub:
#  devtools::install_github("PMassicotte/gtrendsR")
#  library(gtrendsR)

…and we should specify a few key parameters:

# What is the outcome of interest?  
term = "Influenza" #  laboratory-confirmed cases    
# For which country do we want to build the model? 
country_of_interest = "DE" # Germany in ISO_3166-2 (See: https://en.wikipedia.org/wiki/ISO_3166-2)  
# Which language is relevant?
language_of_interest = "de" # German in ISO_639-1 (See: https://en.wikipedia.org/wiki/List_of_ISO_639-1_codes)  
# What the relevant time horizon (i.e. the time span we have data for)
from = as.Date("2010-07-31") # Start
to = as.Date("2017-07-31") # End   
# How do we split the data into training and test data?
split.at = as.Date("2016-08-01") 
# --> Training data:  2010-07-31 - 2016-08-01
# --> Test Data:      2016-08-02 - 2017-07-31

Get data

Outcome data

We can download German influenza incidence data (cases per per 10,000) from the Robert Koch Institute, from August 2010 until August 2017. We uploaded a spreadhseet to github, so it is easier to access. However, you can go to Survstat and customize your data query. The data do not come in the format we need, so we have to re-arrange it a bit.

# Load data from github repository
influenza.de = read.csv("https://raw.githubusercontent.com/projectflutrend/pft.2/master/outcome%20data/RKI.data.de.csv")
names(influenza.de) = c("date","y") 
head(influenza.de)

# Re-formatting 'date', and removing irrelevant data points
influenza.de$date = paste(sub("w","W", influenza.de$date),"-1", sep="")
influenza.de$date  = ISOweek2date(influenza.de$date)
influenza.de = influenza.de[influenza.de$date>as.Date("2010-07-31") & influenza.de$date<=as.Date("2017-07-31"),]
# The RKI data does not report data during summer months, in which there are literally no influenza cases. However, we want to have a model that predicts those zero-months as such, as well. Thus, we will fill the gap and set the cases for the inserted weeks to zero
reference =data.frame(date = seq(from=min(influenza.de$date),
                                 to=max(as.Date(influenza.de$date)),
                                 by=7))
influenza.de = merge(reference,influenza.de,by="date",all=T)
influenza.de$y[is.na(influenza.de$y)] = 0
# Plotting the data
ggplot(influenza.de,aes(x=date,y=y)) + 
  geom_line() + 
  theme_bw() + 
  geom_line() +
  annotate("rect", xmin=split.at, 
           xmax=max(influenza.de$date)+2, 
           ymin=min(influenza.de[,2])-2, 
           ymax=max(influenza.de$y), 
           fill="orange",
           alpha=0.2) +
  annotate("text", 
           x=median(influenza.de$date), 
           y=25, 
           label= paste("Split between training/test data:",split.at), 
           size=4,col="orange") +
  ylab("Incidence per 10,000") +
  ggtitle("Influenza in Germany")

The figure shows the influenza seasons 2010/2011 to 2016/2017, in Germany. More specifically, it shows ‘laboratory-confirmed cases’. It should be noted that this type of outcome is probably not ideal to use as a gold standard for training a digital surveillance system that is based on symptoms: If people feel like having the flu, they will look up influenza on Google and Wikipedia. Flu-like symptoms that occur during the summer months are, however, not caused by the influenza virus- this is why the Robert Koch Institute doesn’t report any data for these months and why we have to assume a flat line at zero. Our prediction model has to deal with the challenge of distinguishing between people who have the real influenza and people who only have influenza-like symptoms.

As a comparison: pattern of influenza-like illness in Germany, as reported by the Robert Koch Institute. Influenza-like illness is probably more appropriate as an outcome for nowcasting models, as it is a better reflection of the incidence of influenza-like symptoms. However, for Germany, this data is not publicly available.

The ‘null-model’

Before we go on to build a prediction model, we should quickly define what we will compare our model against. It would not be fair to compare a prediction model against a flat line at the mean, or at zero cases. This would be like a clinical trial with no comparator. So, let’s compare against placebo: Our model should at least be better than a ‘naiv’ forecast, which only knows what has happend in the past, and extrapolates the pattern into the future (a future that is the present to us). To build this Null-Model, we will use the prophet package that was build by people at Facebook. It is supposed to automatically detect changes in trends, fit a piecewise linear growth curve and add a yearly seasonal component, modeled using Fourier series.

For the present data, we will split the data set into training (seaons 2010/2011-2015/2016) and test data sets (season 2016/2017), we fit a forecast model using the training data and ask for a forecast for the next year. Subsequently, we can compare the predictions with the actual data from the test data set.

# Spliting the data for the forecast model into train and test data set
train.influenza.de = influenza.de[influenza.de$date<split.at,]
test.influenza.de = influenza.de[influenza.de$date>=split.at,]
# 'Propheting'
m <- prophet(df=data.frame(ds = train.influenza.de$date,
                           y=train.influenza.de$y),
             growth = "linear",
             yearly.seasonality = T,
             weekly.seasonality = F)

Initial log joint probability = -6.91857
Optimization terminated normally: 
  Convergence detected: relative gradient magnitude is below tolerance

future <- make_future_dataframe(m, periods = 365)
forecast <- predict(m, future)
p = plot(m, forecast) +
  geom_point(data=test.influenza.de, aes(x=date,y=y),col="pink",lty="dashed", size=2)
p

The prophet-forecast algorithm tries to find a pattern in the available data (2010/2011-2015/2016) and extrapolates it into the future (2016/2017). We can see that the 2016/2017 influenza season was unusually strong, compared to the last six seasons. The forecast model would have underpredicted the number of cases. Furthermore, the season started ealier than usual. So, the forecast model predicted the onset, peak and end of the season a few weeks earlier than they occured. Later, we can use this forecast and compare it to what our Nowcast-model predicts. There are probaly better models to forecast influenza, but it is still better than nothing.

Google Correlate

Google Correlate is a tool that can identify search queries that are highly correlated with any time series data you upload (Even for randomly generated numbers it will find a good match). It was a building block of Google Flu Trend. In addition, you can also identify queries that are correlated with another query: For example, find queries that are related to “Influenza” queries. See this paper for more info on how Google Correlate and Trends works, and this guide on how to use it.

Unfortunately, there is no Google Correlate API available for R, so you need to open the website in your browser and upload your weekly outcome data manually. Google gives you 100 correlated search queries within the country you select (Not all countries are available) and you can download the results as a .csv file. What you need for that is a) A Google Account and b) A spreadsheet with your data in a specific format. Here, we don’t want to leak any information from the test data into the training data- So we will ask for correlated queries only for the trained data set, withholding any patterns seen in the test data set. It’s worth noting that there is a “shift series” button. It shifts your data one week in time, if you think there might be a delay between people using Google and the actual reporting of influenza cases (Only shift your data into the past, the predictors should be collected during the week the cases were reported or BEFORE the cases were reported, not afterwards).

For some reason, Google Correlate won’t provide data for the entire time span. It only gives data points untill 2017-03-12. Therefore, we need to extract only the keywords and download the datapoints from google trends.

# Prepare the outcome data in the right format
g.cor.influenza.training.upload = influenza.de[influenza.de$date<split.at,]  
# Adjusting week format: making Sunday the first day of the week
g.cor.influenza.training.upload$date = g.cor.influenza.training.upload$date-1  
# Saving the file in your default folder
# write.table( g.cor.influenza.training.upload, col.names=FALSE,row.names = FALSE,sep=",",file="/FILL_IN_A_PATH/g.cor.influenza.training.upload.csv")  
### --> Now, go to https://www.google.com/trends/correlate/ 
###     Login, upload the spreadhseet, and download results

Upload data to Google Correlate	Google Correlate Results

# We have put the Google Correlate data for this example on Github:
github.url = "https://raw.githubusercontent.com/projectflutrend/pft.2/master/input%20data/correlate-g_cor_influenza_training_upload.csv"
# The first 10 lines are text, so we want to skip this part
g.cor.results = read.csv(skip=10,github.url) 
# extracting names, except 'date' and 'y'
g.cor.keywords = names(g.cor.results)[-c(1,2)] 
# R inserted "." for spaces, we have to undo this:
g.cor.keywords = gsub("\\."," ",g.cor.keywords)
g.cor.keywords[1:20] # Showing the first 20 keywords

 [1] "influenza a"                  "influenza"                   
 [3] "virusgrippe"                  "j11 1"                       
 [5] "grippe fieber"                "grippe verlauf"              
 [7] "influenza schnelltest"        "influenza inkubationszeit"   
 [9] "echte grippe"                 "grippe bei kindern"          
[11] "grippe husten"                "grippe influenza"            
[13] "ist grippe ansteckend"        "influenza grippe"            
[15] "symptome grippe"              "j11 1 g"                     
[17] "grippe wie lange"             "symptome influenza"          
[19] "wie lange dauert eine grippe" "verlauf grippe"              

If you wonder what “J11 1” might be - It’s the ICD-10 code for Influenza. Doctors print it on the letters for sick leaves. It seems as if some patients were curious what the doctor had diagnosed them with. Other terms refer to [5] fever, [10] Flu in Children, [11] influenza coughing or [17] Influenza how long

Google trends

Google Trends is “a public web facility of Google Inc., based on Google Search, that shows how often a particular search-term is entered relative to the total search-volume across various regions of the world, and in various languages” - (from Wikipedia). It can also be used to identify related search queries (“People who searched for this, also searched for this” - notice the difference between Google Trends and Correlate!)

First, we are going to download the data for the 100 keywords we have identified using Google Correlate. Second, while Google Correlate will most often provide you with good predictors, there might be additional value in taking into account predictors from Google Trends as well, even if only to control for confounding or media-generated over-prediction. Thus, we will download search query statistics for ‘influenza’, and for 15 related queries. Third, we retrive data for ‘influenza’-queries in Google News. If you like, you can expand your requests and increase the numbers of predictors by also retrieving related keywords from related keyword, feed related keywords into Google Correlate, add more News-queries etc.

Cave: Google Correlate only gives you weekly data for time spans < 5 years. As a possible work-around, we split our time series into smaller chunks and download them separately. But because of the way the data is generated, overlapping time series do not match 100%, even after adjusting the scale. Also, it is not clear whether results will be reproducible, as Google might show slightly different figures in new queries.

# Identifying 15 related search queries
  # We have already set these parameters: 
  # term = "influenza"; country_of_interest="DE"; from="2010-07-31"; to="2017-07-31"; status= 1
  google_primer = gtrends(keyword=term,
                          geo=country_of_interest,
                          time=paste(from,to),
                          gprop ="web")
  google_related = google_primer$related_queries$value[google_primer$related_queries$related_queries=="top"]
  g.trends.keywords = c(term,google_related)
  print(g.trends.keywords)

 [1] "Influenza"                      "influenza grippe"              
 [3] "grippe"                         "virus influenza"               
 [5] "rki"                            "influenza rki"                 
 [7] "influenza symptome"             "influenza impfung"             
 [9] "symptome"                       "inkubationszeit"               
[11] "influenza inkubationszeit"      "hämophilus influenza"          
[13] "influenza 2017"                 "influenza 2016"                
[15] "fieber"                         "influenza 2015"                
[17] "influenza 2013"                 "influenza a virus"             
[19] "influenza ansteckend"           "haemophilus influenza"         
[21] "influenza wie lange ansteckend" "grippewelle"                   
[23] "dauer influenza"                "influenza viren"               
[25] "influenza impfstoff"            "was ist influenza"             

The loops to download the search queries are messy and rather long. We will just load them from github, where you can have a closer look at them, if you want.

At the moment, functions only work for the given time span

dim(google.input.data)

[1] 366 128

head(google.input.data[,1:4])

If you want to reproduce this example and avoid long downloading times, you can download the data here

Wikipedia data

McIver & Brownstein made a strong case for the utility of Wikipedia page view data in influenza prediction models. They showed that “Wikipedia usage accurately estimated the week of peak ILI activity 17% more often than Google Flu Trends data and was often more accurate in its measure of ILI intensity”. Combining Google data with Wikipedia may have further advantges.

We can access Wikipedia page view statistics via it’s API (for data from October 2015 until today) and via Wikishark (for data from January 2008 until December 2016)

# On a side note: On Wikipedoa, pages with similar content are linked with each others across different languages. The [Wikipediatrend package](https://cran.r-project.org/web/packages/wikipediatrend/README.html) offers a convenient way to identify these corresponding pages: For 'Influenza', there are about 80 similar pages. While finding the German word for influenza (it's 'influenza'), is not too impressive, automatically finding Wikipedia disease pages in Arabic and Japanese can be quite useful. Wikipedia refers to languages in [ISO 639-1 code](https://en.wikipedia.org/wiki/List_of_ISO_639-1_codes).
  
# The 'Influenza' page in other languges
wikipediatrend::wp_linked_pages( page= "Influenza",lang="en")[1:10,2:3]

   lang    title        
1  af      Griep        
2  als     Influenza    
3  ar      إنفلونزا     
4  an      Gripe        
5  roa-rup Aremi        
6  as      ইনফ্লুৱেঞ্জা    
7  ast     Gripe        
8  gn      Mba'asypararã
9  ay      Jurma_usu    
10 az      Qrip         

To download Wikipedia page view data, we have to turn to two sifferent sources of data. Wikipedia has its own API, which allows fast and convenient open access. However, the API is only able to retrieve data for > October 2015. Statistics for prior dates are stored in huge files, showing page views per hour in every language. Fortunately, there is a private project, called wikishark, which we can quickly use to download relevant data per week (There used to be another server with an API: stats.grok.se, but it appears to be down since July 2017). For the two time periods, Wikipedia has used different metrics to count page views and to detmerine individual page visitors. Thus, there might be some discrepancies. It is also important to note that Wikipedia is not static. New pages are created and old pages may be changed or merged, i.e. for some pages, that do exist today, we won’t find any meaningful page view statistics in 2010 - So, expect (and ignore) some errors when downloading the data. For further information about limitations and potential influences on view counts, see the discussion and related pages on Wikipedia.

First, we download page view statistics for ‘influenza’. From this page, we can also extract all the links, which refer to other Wikipedia pages (E.g. there are links to pages about Aspirin, allergy, bronchitis etc.), and we can also get the names of the pages which link to the influenza page (Treatments for influenza link to the influenza page, for example, but most ‘backlinks’ are more loosely associated). If you like, you can also add Wikipedia pages manually.

The functions to identify linked pages and to download the statistics are rather long and messy. We cut it short and recommend just loading the functions from github, without going through the code. If you like, you can view the code there.

At the moment, functions only work for the given time span

cat("Data for",length(wikipedia.input.data) -1, "Wikipedia pages were downloaded. \n")

Data for 408 Wikipedia pages were downloaded. 

head(wikipedia.input.data[,1:4])

If you want to reproduce this example and avoid long downloading times, you can download the data here

Merging

Now that we have all the data we need, we can merge the files into one dataframe.

# Combining outcome, wikipedia, google trends and google correlate
influenza.de$date = ISOweek(influenza.de$date ) 
# Merging by week (avoiding any Monday/Sunday or other day issues)
df.full = merge(influenza.de,google.input.data, by="date")
df.full = merge(df.full,wikipedia.input.data, by="date")
# Setting date back to a date
df.full$date = ISOweek2date(paste(df.full$date,"-1",sep="")) # 
cat("Full data set:",dim(df.full)[1], "Weeks and",dim(df.full)[2]-2,"Predictors")

Full data set: 346 Weeks and 535 Predictors

This data set can also be downloaded

---

Pre-processing

Before the data can be used to build models, we have to split it into training and test data, deal with missing observations, remove predictors with too little variance, transform skewed predictors, think about multi-colinearity and bring all predictors to a common scale.

Splitting the data

Splitting the data into predictor and outcome, as well as training and testing data sets.

# Remeber: 'split.at' defines the date which separates training and test/evaluation data
split = which(df.full$date<split.at) 
df.train = df.full[split,-c(1,2)] # Predictor training data set
y.train = df.full[split,c(2)] # Outcome for training data set
date.train = df.full[split,c(1)] # Date, not a predictor but useful for plotting
df.test  = df.full[-split,-c(1,2)] # Predictors for testing/evaluation data set
y.test = df.full[-split,c(2)] # Outcome for testing data set
date.test = df.full[-split,c(1)] # date for test data set

Missing data points

Furthermore, we remove predictors that have too many (say, >10%) missing values (see figure below). This method will be applied to the test data set, as well- If we wouldn’t, we could end up with predictors that have too many missing values in the test data set to give a proper prediction. Then, we would need to re-run the whole model again without the missing test-predictor, anyway.

When less than 10% of cases are missing, we are going to impute the missing values for the training and test data set separately, using an exponatially weighted moving average.

# Removing features with >10% NAs
  # in training 
  sum.NA.train = as.numeric(lapply(df.train,function(x){sum(is.na(x))})) 
  sum.NA.train = sum.NA.train > length(df.train[,1]) * 0.1 
  if(sum(sum.NA.train)>0){
  df.train = df.train[-which(sum.NA.train)]
  df.test = df.test[which(colnames(df.test) %in% colnames(df.train))]}
  # and test data separately
  sum.NA.test = as.numeric(lapply(df.test,function(x){sum(is.na(x))}))
  sum.NA.test = sum.NA.test > length(df.test[,1]) * 0.1 
  if(sum(sum.NA.test)>0){
  df.test = df.test[-which(sum.NA.test)]
  df.train = df.train[which(colnames(df.train) %in% colnames(df.test))]}
  
# Imputing remaining NAs
  df.train = na.ma(df.train , k = 3, weighting = "exponential") 
  df.test = na.ma(df.test , k = 3, weighting = "exponential") 

Low variance predictors

Furthermore, we want to remove predictors that have near zero variance. This means they have few unique values. Take for example the Wikipedia page on “Samuel Warren Abott”: The page was created only in 2015, so in the time before, the page had mostly 0 page views (Only rarely someone tries to access a Wikipedia pages that are non-existent). Prediction models can be adversely affected by these near zero variance predictors. You might also consider removing predictors that have a very low activity, because results could be very unstable if models put much weight on them (a random increase in interest could distort our model).

# Example: near zero variance 
nearzero.sample = nearZeroVar(df.train,freqCut = 95/5 , uniqueCut = 33)[11] 
ggplot(df.train,aes(x=date.train ,y=df.train[,nearzero.sample])) +
  geom_line() +
  ylab("activity") +
  ggtitle(paste("Near zero variance: ",names(df.train)[nearzero.sample]))

# Removing features with near zero variance
  # identify near zero-variance predictors [only in df.train!]
  nearZeroVar = nearZeroVar(df.train,freqCut = 95/5 , uniqueCut = 25) 
  if(sum(nearZeroVar)>0){
  df.train = df.train[,-nearZeroVar] 
  df.test = df.test[which(colnames(df.test) %in% colnames(df.train))]}

*It is not clear what is going on with the article apout ’Gebr._Klingenberg’ (a German company). Between 2010 and the end of 2015, it was visited only during four weeks (2,1,1,2 page views). Thus, it shows very little variance with regard to the total time span. After October 2015, with the introduction of the new pageview API, we see a drastical change and much more variation. It is not clear what caused this strange behaviour*

Center - Scaling - Skewness reduction

Predictors should be scaled (each value of a predictor is divided by its standard deviation) and centred (All predictor have a mean of zero). This can improve the performance of some of the modelling techniques and makes it easier to compare predictors’ influence. Be aware, though, that the information, about how many absolute visits a Wikipedia page had, will be lost and that can also be more difficult to interpret the results.

Furthermore, many of our predictors are heavily right-skewed. In order to improve the performance of some of our modelling techniques, we will use caret’s ‘BoxCox’ function which automatically determines the appropriate transformation to make the data ‘more normal’.

# Scaling, centering, transofrmation and imputation of remaining NAs by K-nearest neighbours
  preprocess.df.train = preProcess(df.train, method=c("scale","center","BoxCox"))
  df.train = predict(preprocess.df.train, newdata = df.train)
  df.test = predict(preprocess.df.train,newdata = df.test)

Mutli-colinarity

Many of our predictors are strongly correlated with each other - especially the Google Correlate predictor set (see figure below). To takle this issue, you might consider removing some of the predictors, or use principal component analysis (PCA), to reduce the number of predictors by creating a smaller set of uncorrelated components. However, PCA is unsupervised, i.e. it is not guided by the outcome variable when reducing dimensions (in contrast to PLS, see below). Since we did not observe any improvment by using PCA with the present data set, we won’t cover it here.

Most predictors in the Google Correlate Predictor Set are correlated with each other. “Grippe Verlauf” and “Influenza Inkubationszeit” show a correlation of r=.85 for example, and the mean correlation is around r=.8 over the selected 11 predictors

# If you want to use PCA to tackle this issue in R:
 PCA.proces = preProcess(df.raw, method=c("pca"),  # principal component analysis
                          thresh=0.95)             # a cutoff for the cumulative percent of variance                                                    # to be retained by PCA

 df.pca = predict(PCA.proces, newdata = df.raw)

---

Model building

Overview

We will run several types of models:

Models	Resources
partial least squares
ridge regression
lasso regression	Tibshirani for a technical or this website and video for a more practical explanation
multivariate adaptive regression splines
support vector machine
single trees
boosted trees
bagged trees
random forest
Cubist
Neural Network

Models will be judged by the accuracy of their predictions root-mean-squared-error (rmse), within the cross-validation set. We will evaluate the performance of the model(s) with the lowest rmse within the test data set.

Time series cross validation

Usual split-sample cross validation is problematic in time series data. A random split of the data, as well as splitting the data based on season-year, would leak information from the future to the past. This is a critical issue in predictors that loose or gain predictive power over time. Consider the following simple example: When we train a model with data from 2010 and 2012 and validate it with data from 2011, information about the predictive power of the keyword ‘influenza 2010’ (which is very high in 2010 but low in 2011) is leaked from 2012 (where it is also low) to 2010.

Therefore, one should use time series cross validation. This means, we will train a model using the data from the first data point until week x and validate the results using data from the next k-weeks. In each validation round, we will move k weeks forward, while the origin stays the same. The amount of data the model is being trained on increases with every round (A model is trained using data from 2010-W30 to 2011-W30 and validated with data from 2011-W31 and W32, in the next round, the model is trained from 2010-W30 to 2011-W32 and validated with data from 2011-W33 and W34, and so on). You could also consider using a moving origin (i.e. the amount of training data stays the same), if you think online behaviour from six years ago is not relevant any more and models should only be based on data from the last two seasons, for example. Small k’s provide more accurate rmse results, but it also increases computational time. In this example we chose a k of 13.

# Setting rolling forward cv with fixed origin
controlObject <- trainControl(method = "timeslice",
                              initialWindow = 52,  # First model is trained on 52 weeks (x)
                              horizon = 52,        # Validation weeks (k)
                              fixedWindow = FALSE, # Origin stays the same
                              allowParallel = TRUE)# Paralel computing can speed things up

Parallel computing

Parallel computing can shorten the time it takes to run these computationally intensive models signficantly.

# parallel computing
no_cores <- detectCores() - 1  
cl <- makeCluster(no_cores, type="FORK")
registerDoParallel(cl)  

Modelling

I have only run a few of models, yet. Takes so much time… Note: Save models on github…

start.time = Sys.time()

# partial least square
  M.pls = train(y= y.train ,
              x = df.train,
              method = "pls",
              tuneLength = 20,
              trControl = controlObject)

# ridge regression (enet)  
  # ridge grid
  ridgeGrid <- expand.grid(.lambda = c(0, .001, .01, .1),.fraction = seq(0.05, 1, length = 10))  
  # model
  M.ridge  = train(y= y.train ,
                 x = df.train,
                 method = "enet",
                 tuneGrid = ridgeGrid,
                 trControl = controlObject)

# lasso regression (glmnet)
  # lasso grid
  lassoGrid <- expand.grid(.alpha = c(0,.2, .4, .6, .8, 1),.lambda = seq(.01, 2, length = 20))
  # Model
  M.lasso <- train(y= y.train ,
                     x = df.train,
                     method = "glmnet",
                    family = "poisson", # gaussian?!
                    tuneGrid = lassoGrid,
                    trControl = controlObject)

# multivariate adaptive regression splines (earth)
  # mars grid
  marsGrid <- expand.grid(.degree = 1, .nprune = 2:25)
  # Model
  M.mars=train(y= y.train ,
          x = df.train,
          method = "earth",
          tuneGrid = marsGrid,
          trControl = controlObject)
  
    # varImp(M.mars)

# radial support vector machine (svmRadial)
  # rsvm grid
  rsvmGrid <-  expand.grid(sigma= 2^c(-25, -20, -15,-10, -5, 0), C= 2^c(0:4))
  # Model
  M.rsvm = train(y= y.train ,
              x = df.train,
              method = "svmRadial",
              tuneLength = 25,
              tuneGrid =rsvmGrid, 
              trControl = controlObject)

# Single regression trees 1 (rpart)
  M.rpart = train(y= y.train ,
                x = df.train,
                method = "rpart",
                tuneLength = 20,
                trControl = controlObject)
  
# Single regression trees 2 (cpart)
  M.ctree = train(y= y.train ,
                x = df.train,
                method = "ctree",
                tuneLength = 20,
                trControl = controlObject)

# Bagged trees (treeBag)
  M.bagtree = train(y= y.train ,
                      x = df.train,
                      method = "treebag",
                      trControl = controlObject)
  
# Random forest (rf)
  M.rf = train(y= y.train , ### takes too long!!!
              x = df.train,
              method = "rf",
              tuneLength = 10,
              ntrees = 500,
              importance = TRUE,
              trControl = controlObject)
  
# Boosted tree (gbm)
  # boostree grid
  boostGrid = expand.grid(.interaction.depth = seq(1, 7, by = 2), .n.trees = seq(100, 1000, by = 50),.shrinkage = c(0.01, 0.1),.n.minobsinnode = c(10))
  # Model
  M.boost = train(y= y.train ,
                  x = df.train,
                  method = "gbm",
                  tuneGrid = boostGrid,
                  verbose = FALSE,
                  trControl = controlObject)
  
# Cubist (cubist)
  # cubist grid
  cubistGrid <- expand.grid(.committees = c(1, 25, 50, 75, 100,125,150,175,200),.neighbors=c(0,1,5,9))
  # Model
  M.cubist = train(y= y.train ,
                   x = df.train,
                   method = "cubist",
                   tuneGrid = cubistGrid,
                   trControl = controlObject)
  
# neural network (avNNet)
  # nnet grid
  nnetGrid <- expand.grid(.decay = c(0.001, .01, .1),.size = seq(1, 27, by = 2),.bag = FALSE)
  # model
  M.nnet = train(y= y.train ,
                x = df.train,
                method = "avNNet",
                tuneGrid = nnetGrid,
                preProcess = "pca", # nnet are often affected by multicolinearity
                linout = TRUE,
                trace = FALSE,
                maxit = 500,
                trControl = controlObject)
  
  finished.time = Sys.time() - start.time
  cat(finished.time)

for(i in 1:length(models.de$result.list)){
  tryCatch({
  models.de$eval.list[[i]] = pft_eval_model(models.de$result.list[[i]])
    names(models.de$eval.list)[i] = names(models.de$result.list)[i]},
  error = function(e){cat(t ,": Error \n")})
}

4 : Error 

---

Results

Comparing model results

Evaluating the models individually

How you should read the figure below:

Each row represents one type of model. On the left side, you see the results from the model tuning/cross validation. The rmse refers to the mean error within the validation set. The tuning parameters with the lowest rmse are selected for the final model. On the right side, you see the comparison between predicted and observerd values, of the final model, within the complete training data set and the test data set. We expect to see a lower error within the training data set (right plot, orange), than in the test data set (right plot, red) or the validation set (left side). Furthermore, the mean rmse within the validation sets should be comparable to the rmse within the test data set. The most important

NOT INCLUDED, YET:

                         plot(models.de$result.list[[4]] ),
                         models.de$eval.list[[4]]$plots$pred.plot,
                         plot(models.de$result.list[[5]] ),
                         models.de$eval.list[[5]]$plots$pred.plot,
                         plot(models.de$result.list[[6]] ),
                         models.de$eval.list[[6]]$plots$pred.plot,
                         plot(models.de$result.list[[7]] ),
                         models.de$eval.list[[7]]$plots$pred.plot,

Model selection

Nowcast model versu Forecast model

Influenza-Nowcasting

---

Discussion

Change point analysis

more data sources: twitter, etc etc

transfer the approach to other settings and see what the results are. To the best of our knowledge, digital epidemiology has not been tried in most countries around the world. Results might differ considerably, depending on the degree of digitalization, information seeking behaviour, the sources of information, etc etc.

reporting

Text

…render(“input.Rmd”, html_document(code_download=TRUE))